% Scrip to import CrowdMag data and analize it

%load data
fid = fopen('/Users/manojnair/data/crowdmag/export_May_4_2015.dsv','rt')
[A,B] = textscan(fid,'%s%s%f%f%f%f%f%f%s%s%s%f%f','Delimiter','|','Headerlines',1);

% columns in the dsv file
% DEVICE|OBSTIME|LONGITUDE|LATITUDE|ACCURACY|ELEM_X|ELEM_Y|ELEM_Z|OBJECTID|SHAPE|VALID|ALTITUDE|INTENSITY
%% convert timestamp string to fday

timestampcell = cell2mat([A{2}]);
timestampcell(:,19:28) = []; % remome decimal seconds
timestampcell(1,:)
fday = datenum(timestampcell,'dd-mmm-yy HH.MM.SS PM');

%% get devices matrix

device = [A{1}];

%% get the data
c_latitude = [A{4}];
c_longitude = [A{3}];
c_location_accuracy = [A{5}];
c_h = [A{7}];
c_z = [A{8}];


%%
ncell = 1;
delta_lat = 0.25;
longitudearray = -179:delta_lat:179;
latitudearray = -80:delta_lat:80;

wmm_ngdc = '/Users/manojnair/projects/wmm2015_validation/WMM2015_FINAL.txt';


for j = 1:length(latitudearray),
    tic;
    for i = 1:length(longitudearray),
        
        longitude = longitudearray(i);
        latitude = latitudearray(j);
        
        L = c_latitude >= latitude & c_latitude <= latitude + delta_lat ...
            & c_longitude >= longitude & c_longitude <= longitude + delta_lat;
        
        
        if any(L),
            
            % get WMM2015 values
            
            geoc_lat = geod2geoc(latitude + delta_lat/2, 0, 'WGS84'); % geocentric lat
            geoc_r = ( geocradius(geoc_lat,'WGS84') - 6371.2*1000 )/1000; % geocentric alt in km
            wmm_2015 = magsynth(geoc_r,geoc_lat,longitude + delta_lat/2 ,2015.0,wmm_ngdc );
            
            
            data_cell(ncell).wmmx = wmm_2015(1);
            data_cell(ncell).wmmy = wmm_2015(2);
            data_cell(ncell).wmmz = wmm_2015(3);
            data_cell(ncell).z = median(z(L));
            data_cell(ncell).ndata = sum(L);
            data_cell(ncell).h = median(h(L));
            data_cell(ncell).hstd = std(h(L));
            data_cell(ncell).zstd = std(z(L));
            data_cell(ncell).lat = latitude + delta_lat/2;
            data_cell(ncell).lon = longitude + delta_lat/2;
            data_cell(ncell).devices = unique(device(L));
            data_cell(ncell).ndevices = length(unique(device(L)));
            
            ncell = ncell + 1;
            
            %fprintf('%5.2f %d \n',longitude);
        end;
        
        
    end;
    
    fprintf('%5.2f %s \n',latitude, toc);
    %break;
end;


%% Get WMM on a global grid


wmm_ngdc = '/Users/manojnair/projects/wmm2015_validation/WMM2015_FINAL.txt';

ilat = 1;
for lat = -89.5:89.5,
    ilon = 1;
    
    for lon = -179.5:179.5,
        
        geoc_lat = geod2geoc(lat, 0, 'WGS84'); % geocentric lat
        geoc_r = ( geocradius(geoc_lat,'WGS84') - 6371.2*1000 )/1000; % geocentric alt in km
        wmm_2015 = magsynth(geoc_r,geoc_lat,lon,2015.0,wmm_ngdc );
        
        WMM_X(ilat,ilon) = wmm_2015(1);
        WMM_Y(ilat,ilon) = wmm_2015(2);
        WMM_Z(ilat,ilon) = wmm_2015(3);
        
        ilon = ilon +1;
        
    end;
    ilat = ilat + 1;
    fprintf('%d\n', lat);
end;

%% Plot
ndata = [data_cell.ndata];
latitude = [data_cell.lat];
longitude = [data_cell.lon];
z = [data_cell.z];
wmmz = [data_cell.wmmz];
wmmh = sqrt([data_cell.wmmx].^2 + [data_cell.wmmy].^2);
h = [data_cell.h];
WMM_F = sqrt(WMM_X.^2+WMM_Y.^2+WMM_Z.^2);;
f =sqrt(z.^2+h.^2);

L = ndata < 20;

worldmap('world')
load coast
plotm(lat,long);
scatterm(latitude(~L),longitude(~L), ones([1,sum(~L)])*200, z(~L),'filled');
caxis([-65000 65000]);
colorbar


L = ndata < 0;


worldmap('world')
load coast
plotm(lat,long);
scatterm(latitude(~L),longitude(~L), ones([1,sum(~L)])*100, (ndata(~L)),'filled');
%caxis([-65000 65000]);
colorbar


%% plot Z from crowdMag and WMM

L = ndata < 1000 | ( abs(latitude) <1 & abs(longitude <1))  | ndevices < 4;

worldmap('world')
load coast
plotm(lat,long);
scatterm(latitude(~L),longitude(~L), ones([1,sum(~L)])*500, z(~L),'filled');
[c,h] = contourm(WMM_Z, [1, 89.5, -179.5],[-650000:10000:65000]);
clabelm(c,h);
caxis([-65000 65000]);
colormap(jet);
h = colorbar;
set(h,'FontSize',16);


L = ndata < 1000 | ( abs(latitude) <1 & abs(longitude <1))  | ndevices < 4;

worldmap('world')
load coast
plotm(lat,long);
scatterm(latitude(~L),longitude(~L), ones([1,sum(~L)])*500, f(~L),'filled');
[c,h] = contourm(WMM_F, [1, 89.5, -179.5],[0:10000:75000]);
clabelm(c,h);
caxis([0 75000]);
colormap(jet);
h = colorbar;
set(h,'FontSize',16);


L = ndata < 1000 | ( abs(latitude) <1 & abs(longitude <1))  ;
zstd = [data_cell.zstd];
worldmap('world')
load coast
plotm(lat,long);
scatterm(latitude(~L),longitude(~L), ones([1,sum(~L)])*500, log10(zstd(~L)),'filled');
%clabelm(c,h);
%caxis([0 75000]);
colormap(jet);
h = colorbar;
set(h,'FontSize',16);
%% Temporal changes

% get the Nexus 4 data
latitude = 39.55;
longitude = -105.0;
delta_lat = 0.2;

L = strcmp(b(:,1),'Nexus 4') &  a(:,2) >= latitude & a(:,2) <= latitude + delta_lat ...
    & a(:,1) >= longitude & a(:,1) <= longitude + delta_lat;

% Get F data from this station

fday_nexus = fday(L);

f_nexus = sqrt(a(L,6).^2 + a(L,5).^2);
h_nexus = a(L,5);
z_nexus = a(L,6);

%sort time series
[y,i] = sort(fday_nexus);

%compared with Boulder F (12/3-12/8). Nexus variation is about 10 times as large as
%Boulder

%make histogram

%Filtering is justified since the plot shows data for F
%below this value is very strange
L = f_nexus > 49800;

f_nexus_filtered = f_nexus(L);
h_nexus_filtered = h_nexus(L);
z_nexus_filtered = z_nexus(L);
fday_nexus_filtered = fday_nexus(L);
[y,i] = sort(fday_nexus_filtered);


% iPhone5,2
latitude = 39.5;
longitude = -105.5;
delta_lat = 1;

L = strcmp(b(:,1),'iPhone5,2') &  a(:,2) >= latitude & a(:,2) <= latitude + delta_lat ...
    & a(:,1) >= longitude & a(:,1) <= longitude + delta_lat;


fday_iphone = fday(L);

f_iphone = sqrt(a(L,6).^2 + a(L,5).^2);
h_iphone = a(L,5);
z_iphone = a(L,6);

% sort time series
[y,i] = sort(fday_iphone);

%% Look at any spatial patterns

ndata = [data_cell.ndata];
latitude_all = [data_cell.lat];
longitude_all = [data_cell.lon];
z = [data_cell.z];
wmmz = [data_cell.wmmz];
wmmh = sqrt([data_cell.wmmx].^2 + [data_cell.wmmy].^2);
h = [data_cell.h];


% select grids whith atleast 10,000 data points

L = ndata>10000;

selected_lat = latitude_all(L);
selected_lon = longitude_all(L);
delta_lat = 0.25;
%%
latitude = selected_lat(6) - delta_lat/2;
longitude = selected_lon(6) - delta_lat/2;

% MIT location of interest: axis([ -71.0899  -71.0891   42.3599   42.3610])
% 6, ndata>10000 detla_lat  0.25

L = a(:,2) >= latitude & a(:,2) <= latitude + delta_lat ...
    & a(:,1) >= longitude & a(:,1) <= longitude + delta_lat & a(:,3) < 10;

figure(1);
scatter(a(L,1),a(L,2), ones([1,sum(L)])*100, (a(L,6)),'filled');
caxis([median(a(L,6))-std(a(L,6))*3, median(a(L,6))+std(a(L,6)*3)]);
title('Z');
colorbar;
figure(2);
scatter(a(L,1),a(L,2), ones([1,sum(L)])*100, (a(L,5)),'filled');
caxis([median(a(L,5))-std(a(L,5))*3, median(a(L,5))+std(a(L,5)*3)]);
title('H');
colorbar;
figure(3);
colorbar;
scatter(a(L,1),a(L,2), ones([1,sum(L)])*100, (sqrt((a(L,5).^2 + a(L,6).^2))),'filled');
title('F');
colorbar;

%% find out the phone, data and time for MIT anomaly

L = a(:,2) >= 42.3603 & a(:,2) <= 42.3604  ...
    & a(:,1) >= -71.0898 & a(:,1) <= -71.0896 & a(:,3) < 100;

scatter(a(L,1),a(L,2), ones([1,sum(L)])*100, (a(L,5)),'filled');
caxis([median(a(L,5))-std(a(L,5))*3, median(a(L,5))+std(a(L,5)*3)]);
title('H');
text(a(L,1),a(L,2),b(L,1))
scatter(a(L,1),a(L,2), ones([1,sum(L)])*100, (a(L,5)),'filled');
caxis([median(a(L,5))-std(a(L,5))*3, median(a(L,5))+std(a(L,5)*3)]);
title('H');
text(a(L,1),a(L,2),b(L,2))
scatter(a(L,1),a(L,2), ones([1,sum(L)])*100, (a(L,5)),'filled');
caxis([median(a(L,5))-std(a(L,5))*3, median(a(L,5))+std(a(L,5)*3)]);
title('H');
text(a(L,1),a(L,2),num2str(a(L,5)))
